Ion-ion correlations: an improved one-component plasma correction 
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Based on a Debye-Hiickel approach to the one-component plasma we propose a new free energy 
for incorporating ionic correlations into Poisson-Boltzmann like theories. Its derivation employs the 
exclusion of the charged background in the vicinity of the central ion, thereby yielding a thermo- 
dynamically stable free energy density, applicable within a local density approximation. This is an 
improvement over the existing Debye-Hiickel plus hole theory, which in this situation suffers from 
a "structuring catastrophe" . For the simple example of a strongly charged stiff rod surrounded by 
its counterions we demonstrate that the Poisson-Boltzmann free energy functional augmented by 
our new correction accounts for the correlations present in this system when compared to molecular 
dynamics simulations. 
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The classical one-component plasma (OCP) is an idealized model, in which a single species of ions moves 

C^ homogeneous neutralizing background of opposite charge and interacts only via a repulsive Coulomb potentiall 

£h ' Apart from its applications in plasma physicajD it is also commonly used in soft matter physics &$, one of the 

i simplest passible approaches for modeling correlations when studying polyelectrolytes, charged planesQ'u or charged 

' O • colloidsErEij. The general idea is the following: Compute the OCP free energy as a function of bulk density Ub and use 

" ' this expression in the spirit of a local density approximation (LDA) as a correlation correction for the inhomogeneous 
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system (i.e., ne — ► n(r)). The total excess free energy is the volume integral over the free energy density aud thus 
becomes a functional of n(r). Many alternative and more sophisticated methods based on integral equationslij have 
been developed for treating this correlation problem. Even though they offer results which are in good agreement 
with Monte-Carlo simulations, they do not provide any intuitive insight into the physics governing ionic solutions. 

There is, however- a fundamental problem with the local density approaches: The OCP free energy is not a convex 

function of densityEj. This implies that it cannot be used in a thermodynamically stable way within LDA, since the 

r^\ ' system can lower its total free energy by developing locaL-Lnhomogeneities and increasing its density in one region at 

f^ . the expense of another (disregarding any surface effects)Ej. Once started, this continues as a runaway process and 

7— I ' the overall system collapses to a point. This feature is already seen on the level of the Debye-Hiickel plus hole (DHH) 

OS , approximation^, which is an extension of the original Debye-Hiickel (DH) theoryEEl for tha-s»ecial case of the OCP, 

and the instability it gives rise to has been termed "structuring catastrophe" in this contextOEJ. The proper way for 

avoiding this difficulty thus requires modifications of the one-component plasma model itself. The new theory, referred 

to as the Debye-Hiickel-Hole-Cavity (DHHC) approach, remains simple and can be used within LDA to account for 

correlation effects present in more complex ionic solutions, as will be shown in an example at the end of the paper, 

i-rt . where we compare its predictions to simulational results of a model system. 

Since the necessary changes to DHH will turn out to be surprisingly tiny, it is worthwhile to briefly recall the 

O way in which DHH theory arrives at a free energy for the OCP. For definiteness, we assume a system of N identical 

~J ' point-particles of valence v and (positive) unit charge q inside a volume V with a uniform neutralizing background of 

^* . density vn-g, and dielectric constant e. According to the DH approach, the potential 4> created by a central ion (i.e., 

l>( ' fixed at the origin) and all its surrounding ions results from solving the spherically symmetric Poisson equation 

5— i ' 2 Ait 

d- V 2 </>(r) = <f>"(r) + -<f/(r) = p(r) (1) 

r e 

under the requirement that the charge density is p(r) = qvS(r) at the central ion and that the rest of the mo- 
bile ions rearrange themselves in the uniform background in accordance with the Boltzmann distribution p(r) = 
vqnB exp [— (3vq4>(r)] — vqnji- Combining this with eq. (|lj) yields the nonlinear Poisson-Boltzmann (PB) equation, 
while linearization of the exponential function in the mobile ion density gives p(r) = —En 2 (f)(r)/4ir together with 
the famous Debye-Hiickel solution for the potential, 4>(r) ~ e _Kr /r, illustrating the rearrangement of the other ions 
around the central one in order to screen the Coulomb interaction. Here, k = y/AirinB is the inverse screening length, 
i — £bv 2 , with £b = f3q 2 /s being the Bjerrum length, and j3 = l/fceT. 

The problem with the DH theory is that the condition for linearization is obviously not satisfied for small r, where 
the potential is large — indeed, the particle density becomes negative and finally diverges at the origin. This defect 



was overcome by the DHH theorytH, which artificially postulates a correlation hole of radius h around the central ion 
where no other ions are allowed. In this case the charge density is given by 

n(r) - J qV ^ r) ~ " B ) : r ~ h (1\ 

P[r) ~ \ -en 2 (t>{r)/AiT : r > h. [Z> 

The solution of the linearized PB equation with the appropriate boundary conditions (continuity of electric field and 
potential) yields the potential for both regions in dependence of h, which has to be fixed on physical grounds: At low 
temperatures the electrostatic repulsion dominates and the minimum ion separation essentially becomes the mean 
separation, so ft, = (47T71b/3) -1 / 3 . At high temperatures, the hole size can be estimated by balancing Coulombic and 
thermal energy, which gives h — £. A systematic way to interpolate between these two limits results from excluding 
particles from a region where their potential energy is larger than some threshold. A natural choice is the thermal 
energy k#T, which leads to 

nh = w-l with uj = (1 + 3£k) 1/3 . (3) 

Incidentally, this assumption also gives a continuous charge density across the hole boundary. 

Once the potential at the position of the central ion is knowiL the electrostatic contribution to the Hclmholtz free 
energy density can be obtained by the Debye charging procesala, as was done previously by Penfold et alHi: 
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The presented simple DHH analysis of the one-component plasma theory offers considerable insight into ionic 
systemSr-apd is in good agreement with Monte-Carlo simulationsEi3 when fluctuations on the charge density are not 
relevantclJ. In principle one can attempt to include such fluctuations by applying the bulk density-functional theory 
in a local way, i.e., tt-b - ► n(r). The basic idea is to obtain the density distribution via functional minimization of the 
Helmholtz free energy 



pF OC p[n(r)} = Jd 3 r |n(r) In (n(r)V p ) + /3/dhh [n(r)]} 



(5) 



under the constraint of global charge neutrality (V p represents the volume of a particle.) Yet, this variational process 
does not lead to a well defined density profile, since /dhh(^) asymptotically behaves like — n 4 / 3 at high densities and 
therefore is not a convex function - with the implications already mentioned in the introduction. At small densities, 
however, the free energy density is convex and changes to a concave form only beyond a critical density n* « 7.8618/^ 3 
(see fig. ||). Hence, if during the process of actually computing n(r) such a density is never met, the theory does 
not "realize" its asymptotic instability and gives a finite (yet, meta-s,tabie) answer. It has in fact been applied to 
account for correlations in the case of systems with low ionic strengtlJD&tLil). Assuming the case of aqueous solutions 
(£b = 7.14A) and monovalent ions we find a critical density n* ss 36 mol/1, which clearly is high enough to prevent a 
runaway process to set in. However, already for divalent and trivalent ions we find n* w 0.56 mol/1 and 0.049 mol/1, 
respectively, which are sufficiently low to be realized and thus to trigger a collapse. Notice the strong dependence of 
n* on valence, namely, on the sixth power. 

To circumvent the instabilities occurring at high densitieSnCHroneously attributed to the local density approach itself, 
a number of nonlocal free energies have been prnpnspHll m 1 4P A In these weighted density approximations (WDA) the 
local density is replaced by a spatially averaged quantity. The main problem with these methods is that the choice 
of the weighting function is somewhat arbitrary. In most cases it is obtained by relating the second variation of 
the free energy with the direct correlation function. At this point the WDA requires prior information about this 
function, which is not yet available and thus has to be calculated using different approaches (like, e.g., integral equation 
theories). Whatever choice one takes, it is still (i) quite arbitrary and (ii) leads to a series of approximations which 
(Hi) instead of clarifying the physics tend to obscure it. 

The instabilities present in the local DHH approach can be properly overcome by recognizing that the failure of this 
model is due to the (too strong) requirement of local charge neutrality imposed by the LDA: A local fluctuation leading 
to an increase of particle density implies a corresponding increase in background density; therefore the fluctuation is 
not suppressed by an increase in repulsive Coulomb interactions but quite on the contrary favored by its decrease. A 
natural solution for that problem is to decouple the particle density from the background density and to apply the 
LDA just to the former one. This, however, leads to nonlinearities in the solution of the differential equation which 
spoil the simplicity of the DH and DHH approximations. The most simple solution is to exclude the neutralizing 
background from a cavity of radius a placed around the central ion only, which is already sufficient to control the 
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FIG. 1. Free energy density of the DHH (dashed, eq. (M)) and DHHC (solid, eq. did)) theory as a function of density for 
Bjerrum length Is = 7.14A and monovalent/divalent ions (left/right). The arrows mark the points at which the DHH free 
energy density changes from convex to concave. A particle volume of V p = (5A) 3 was assumed. 



unphysical divergence of the particle density. Even though it does not accounts for excluded-volume effectsE3n3, the 
parameter a can in principle be identified with the diameter of the particles. In addition, the exclusion hole for 
a < r < h is retained in order to account for the electrostatic repulsion between two ions. Consequently, the charge 
density, which for the usual DHH theory is given by eq. (0), has now three regions: 
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(6) 



The solution of the linearized PB equation with appropriated boundary conditions gives the potential in those regions: 
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with the abbreviation 
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(8) 



In order to obtain the old theory in the limit a — > we choose the hole size h to yield the same screening (i.e., the 
same amount of charge within h) as the DHH theory, which results in 



nh = [(u - l) 3 + (ko) 



31 1/3 



(9) 



This expression has four important physical limits: zero/infinite temperature and low/high density. At low tem- 
perature the exclusion hole has maximum size and, like in the DHH case, behaves as h = (3/47T7Jb + a 3 ) 1 / 3 . As 
the temperature is increased, the hole size shrinks, but contrary to DHH theory it does not vanishes and fc->aas 
T — ► oo. At small densities, entropic effects compete with the Coulombic repulsion and h = £ + a; for high densities, 
the exclusion hole decreases but is again limited below and h — > a. Using this prescription for h, the Helmholtz free 
energy can be obtained by Debye-charging the fluid: 
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and where u> is the same as in eq. (||). The integral can be solved numerically for given values of £&, v and a. As in 
the DHH approach fluctuations are taken into account by allowing the density to become local; thus, n(r) is obtained 
by minimizing the free energy from eq. (||) with /dhh replaced by /dhhc as given by eq. (|H]). But differently from 
the DHH theory, the Debye-Hiickcl-Holc-Cavity free energy is a convex function of density and thus applicable within 
a local density approximation. This situation is depicted in fig. R1 where we plotted the previous expression of the 
free energy of the DHH theory together with the improved expression of the DHHC approach. Recall that the DHH 
free energy has a point of inflection at a critical density n* « 7.8618/^ 3 , which makes it unstable at high densities - 
particularly for multivalent ionic correlations, as is demonstrated in the right part of fig. |l|. 

As an example, we apply this free energy as a correlation correction in the theoretical description of the screening 
of a charged rod, which is a simple model of biologically relevant stiff polyelectrolytes like DNA, actin filaments 
or microtubules. Much of the thermodynamic behavior of these molecules is determined by the distribution of the 
counterions around the polyion. As a model system we take a rod of radius r$ and line charge density A = 0.959 q/r® 
embedded in a cell of outer radius R = 123.8 r$ and the complementary values ^b/^o — 3, v — 1 and Ib/tq = 1, f = 3 
have been investigated,.-which on the plain PB level both give a fraction of condensed counterion (in the Manning 
sense) of roughly 65%On3. This system is thus strongly charged and one ejteects ionic correlations to become relevant. 
Indeed, the comparison between the distributions obtained by simulationcj and the ones from PB theory shows that 
the mean-field approach fails in the limit of high ionic strength. In reality the ions do not just interact with the 
average electrostatic field but if an ion is present in a position r, it tends to push away other ions from that point. 
This effect becomes important at high densities, low temperatures and for multivalent ions. 

As discussed above, a simple way to improve PB theory is to extend the density functional to include a term of the 
form ( |l0|) which accounts for the correlations. The configurational free energy for the screened macroion solution can 
be partitioned into two terms: 

F P [n(r)} = F PB [n(r)}+ f^rfBHHcHr)}. (12) 

The first part 

F PB [n(r)} = / d 3 r j k B Tn(r) In (n(r )V P ) + - qv n(r)4>[n(r)] } (13) 

contains the ideal gas contribution of the small ions, the interaction with the macroion potential and the mean-field 
interaction between the counterions. Minimization of this expression under the constraint of global charge neutrality 
gives - together with the Poisson equation - the Poisson-Boltzmann equation. The inter-particle correlations are now 
approximately accounted for by adding an excess free energy, which is the second term in eq. (Q~2h — the DHHC free 



energy in local density approximation. The equilibrium ion distribution minimizing the functional (h2J) is most easily 



s 



found by means of a Monte-Carlo solver, as has been proposed elsewherecJ. The fraction of ions within a distance r, 

i r 

P(r) — — / dr 2nr vqn(r), (14) 

A Jro 

obtained following this procedure is illustrated in fig. |2j. Compared to the plain PB result the simulation shows a 
stronger condensation of ions in the vicinity of the rod, an effect which is more pronounced in the trivalent system. In 
both cases the increased condensation is reproduced by the correlation corrected PB functional from eq. ( |l2| ) . While 
in the case ^bAo — 3, V — 1 the theoretical prediction practically overlaps the simulation, it somewhat overestimates 
correlations in the complementary case Ib/tq — 1, v = 3. It must, however, be noted that the ions in the simulation 
also interacted via a repulsive Lennard- Jones potential, giving them a diameter of roughly tq. The expected reduction 
of particle density resulting from the additional hard core is not accounted for in the presented theory, but could 
easily be included along the lines of Refs.EZl. 

In conclusion, we have shown that the failure of the local density approximation for the one-component plasma is 
due to the asymptotically concave free energy employed by the DHH theory. To eliminate this problem, we introduced 
a DHHC approach in which the uniform background is absent in the immediate vicinity of the central ion, which leads 
to a convex, thermodynamically stable free energy. Moreover, the local density functional theory derived from this 
assumption is able to correctly account for the correlations between small ions in the presence of a strongly charged 
macroion. This was demonstrated for the case of a stiff rodlike polyelectrolyte by comparing the integrated charge 
density to simulation results of the same model. A more detailed investigation of the applicability of the LDA to 
rodlike polyelectrolytes and to charged colloids is postponed to future work. 







10 20 

r/r 



50 100 



10 20 

7 r o 



50 100 



FIG. 2. Counterion distribution function P(r) from eq. (114) for two cylindrical cell models with R/ro = 123.8, A = 0.959 q/ro 
and the values for Bjerrum length and valence as indicated in the plots. The solid line is the result of a molecular dynamics 
simulation while the dotted line is the prediction from Poisson-Boltzmann theory. The increased counterion condensation 
visible in the simulation is accurately captured by the extended PB theory (dashed line) using the correction from eq. (flOl). 
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